Fortran For Fun之python调用fortran

发表于 | 分类于

以球谐函数为例，说明如何使用python调用fortran的接口。

fortran 模块

建立fortran模块sh_value, 其中主要包括连带勒让得多项式plm(l,m,x)函数，以及球谐函数ylm(l,m,theta,phi)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
module sh_value
  implicit none
  public :: Ylm,Plm,alm
  private
  real,parameter :: pi_ = 3.141592653589793
contains
  !===========================================================================
  real pure function Ylm(l,m,theta,phi)
    !< Calculate the value of real spherical harmonics
    integer,intent(in) :: l      !< Degree of Ylm
    integer,intent(in) :: m      !Order of Ylm
    real, intent(in)   :: theta  !<the latitudinal coordinates
    real, intent(in)   :: phi    !the longitudinal coordinates
    real :: coef
    coef = alm(l,m)
    if( m>=0 ) then
       Ylm=coef*Plm(l,m,cos(theta))*cos(m*phi)
       return
    else
       Ylm=coef*Plm(l,-m,cos(theta))*sin(-m*phi)
    end if
  end function Ylm
  !===========================================================================
  real pure function alm(l,m)
    !< coefficient of real spherical harmonics
    integer,intent(in) :: l,m
    integer :: i
    real :: fact
    fact = 1.0
    if(m == 0) then
       alm = sqrt((2.0*l+1.0)/(4.0*pi_))
    else
       do i = l-abs(m)+1,l+abs(m)
          fact = fact*i
       end do
       alm = sqrt((2.0*l+1.0)/(2.0*pi_*fact))
    endif
  end function alm
  !===========================================================================
  real pure function Plm (l,m,x)
    !< Calculate the value of Associated Legendre Polynomial by recursion relation
    integer,intent(in) :: l, m
    real, intent(in) :: x
    real :: Pmm, Pmp1m, fact, x2
    integer :: i
    Pmm = 1.0
    if(m>=0) then
       fact = 1.0
       x2 = sqrt(1.0-x*x)
       do i=1,m
          Pmm = Pmm*(-fact)*x2
          fact=fact + 2.0
       end do
    endif
    Pmp1m=x*(2*m+1)*Pmm
    if(l==m)then
       Plm=Pmm
    elseif(l==m+1)then
       Plm=Pmp1m
    elseif( l>=m+2 ) then
       do i=m+2,l
          Plm=((2.0*i-1.0)*x*Pmp1m-(i+m-1.0)*Pmm)/(i-m)
          Pmm=Pmp1m
          Pmp1m=plm
       end do
    end if
  end function Plm
end module sh_value

使用f2py生成动态库

通过f2py生成动态库 sh.so,然后直接在pyhton中使用, 命令如下

1
f2py -c sh_value.f90 -m sh

勒让得多项式

建立python脚本 plot_pl.py 画出勒让得多项式， python需要先安装 matplotlib 和 numpy 库。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import matplotlib.pyplot as plt
import numpy as np
import sh
lmax = 5
theta = np.arange(0,np.pi,np.pi/100)
x = np.cos(theta)
nn = np.size(x);
sty = ['-r','-c','-g','-m','-b','-k']
plt.figure(figsize=(12, 9))
for l in np.arange(0,lmax+1):
    pl = np.zeros([nn],np.float64)
    for i in np.arange(nn):
        xi = x[i]
        pl[i] = sh.sh_value.plm(l,0,xi)
    plt.plot(x,pl,sty[l],label='p'+str(l))
plt.title('Legendre Polynomial')
plt.xlabel('x')
plt.ylabel('pl(x)')
plt.legend(loc='best')
plt.savefig('../img/legendre_polynomial.png')
plt.show()

运行该程序，得到勒让得多项式图像
pl

连带勒让得多项式

plot_plm.py 画出连带勒让得多项式， python需要先安装 matplotlib 和 numpy 库。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import matplotlib.pyplot as plt
import numpy as np
import sh
lmax = 4
theta = np.arange(0,np.pi,np.pi/100)
x = np.cos(theta)
nn = np.size(x);
sty = ['-r','-c','--c','-g','--g','-.g','-m','--m','-.m','-..m','-b','--b','-.b','-..b','-*b']
plt.figure(figsize=(16, 12))
k = 0
for l in np.arange(0,lmax+1):
    for m in np.arange(0,l+1):
        plm = np.zeros([nn],np.float64)
        for i in np.arange(nn):
            xi = x[i]
            plm[i] = sh.sh_value.plm(l,m,xi)*sh.sh_value.alm(l,m)*np.sqrt(4.0*np.pi/(2.0*l+1.0))
        plt.plot(x,plm,sty[k],label='p'+str(l)+str(m))
        k = k+1
plt.title('Associate Legendre polynomial(normalized)')
plt.xlabel('x')
plt.ylabel('plm(x)')
plt.legend(bbox_to_anchor=(1.01,0.89),loc='best')
plt.savefig('../img/associate_legendre_polynomial.png')
plt.show()

运行该程序，得到连带勒让得多项式图像
pl

球谐函数

plot_ylm.py 画出球谐函数， python需要先安装 mayavi 和 numpy 库

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np
import mayavi.mlab as mlab
import sh
pn = 5
theta, phi = np.mgrid[0:np.pi:101j, 0:2*np.pi:101j]
x = np.sin(theta)*np.cos(phi)
y = np.sin(theta)*np.sin(phi)
z = np.cos(theta)
mlab.figure(size=(1200, 900))
for n in np.arange(pn+1):
    for m in np.arange(-n,n+1):
        s = np.zeros([101,101],np.float64)
        for i in np.arange(101):
            t = theta[i][0]
            for j in np.arange(101):
                p = phi[0][j]
                s[i][j]= sh.sh_value.ylm(n,m,t,p)
        mlab.mesh(abs(s)*x-1.5*n, abs(s)*y+1.5*m, abs(s)*z,scalars=s, colormap='Spectral')
mlab.view(azimuth=45.0, elevation=45.0, distance=30.0)
mlab.savefig('../img/spherical_harmonics.png')
mlab.show()

运行该程序，得到球谐函数图像
pl